#!/usr/bin/perl

use strict;
use warnings;

# ------------------------------------------------------------------------------
# --- Convert BSNP file to FASTA of just the BPs in BED file intervals
# ------------------------------------------------------------------------------

my $usage = "Usage: $0 <in.bed> <in.bsnp.snp.out>\n";

my $in_bed = shift;
chomp $in_bed;

my $bsnp_file = shift;
chomp $bsnp_file;

open (BED, "<$in_bed")
    or die "ERROR: Could not open input BED file, $in_bed.\n";

open (BSNP, "<$bsnp_file")
    or die "ERROR: Could not open input BSNP file, $bsnp_file.\n";

my $reduced_bsnp_fasta = $bsnp_file . ".reduced.fa";
open (FASTA, ">$reduced_bsnp_fasta")
    or die "ERROR: Could not open output BSNP file, $reduced_bsnp_fasta.\n";

BED_LOOP: foreach my $bed (<BED>) {

    print STDERR "Processing BED interval: $bed";
    
    my @bed_info = split /\t/, $bed;

    my $bed_chrom = $bed_info[0];
    my $bed_start = $bed_info[1];
    my $bed_end   = $bed_info[2];
    chomp $bed_end;

    print FASTA '>' . $bed_chrom . ':' . $bed_start . '-' . $bed_end . "\n";
    
    BSNP_LOOP: while (my $bsnp = <BSNP>) {
    
        my @bsnp_info = split /\t/, $bsnp;

        if ($bsnp_info[0] ne $bed_chrom) {
            # Close and reopen BSNP file if BED intervals 
            # are still waiting to be processed
            if (eof BSNP) {
                print STDERR "Reopening BSNP\n";
                close BSNP;
                open (BSNP, "<$bsnp_file")
                    or die "ERROR: Could not open input BSNP file, $bsnp.\n";
            }
            next BSNP_LOOP;
        } elsif ($bsnp_info[1] < $bed_start) {
            next BSNP_LOOP;
        } elsif ($bsnp_info[1] == $bed_end) {
            print FASTA $bsnp_info[4];
            print FASTA "\n";
            next BED_LOOP;
        } else {
            print FASTA $bsnp_info[4];
            next BSNP_LOOP;
        }
        
    }

}

close FASTA;
close BSNP;
close BED;

exit;